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ABSTRACT 

We investigate mode coupling in a two dimensional compressible disc with radial stratifica- 
tion and differential rotation. We employ the global radial scaling of linear perturbations and 
study the linear modes in the local shearing sheet approximation. We employ a three-mode 
formalism and study the vorticity (W), entropy (S) and compressional (P) modes and their 
coupling properties. The system exhibits asymmetric three-mode coupling: these include mu- 
tual coupling of S and P-modes, S and W-modes, and asymmetric coupling between the W 
and P-modes. P-mode perturbations are able to generate potential vorticity through indirect 
three-mode coupling. This process indicates that compressional perturbations can lead to the 
development of vortical structures and influence the dynamics of radially stratified hydrody- 
namic accretion and protoplanetary discs. 
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1 INTRODUCTION 

The recent increased interest in the analysis of hydrodynamic disc 
flows is motivated, on one hand, by the study of turbulent pro- 
cesses, and, on the other, by the investigation of regular structure 
formation in protoplanetary discs. Indeed, many astrophysical discs 
are thought to be neutral or having ionization rates too low to ef- 
fectively couple with magnetic field. Among these are cool and 
dense areas of protoplanetary discs, discs around young stars, X- 
ray transient and dwarf nova systems in quiescence (see e.g. Gam- 
mie and Menou 1998, Sano et al. 2000, Fromang, Terquem and 
Balbus 2002). Observational data shows that astrophysical discs 
often exhibit radial gradients of thermodynamic variables (see e.g. 
Sandin et al. 2008, Issela et al. 2007). To what extent these inho- 
mogeneities affect the processes occurring in the disc is still subject 
open to investigations. It has been found that strong local entropy 
gradients in the radial direction may drive the Rossby wave insta- 
bility (Lovelace et al. 1999, Li et al. 2000) that transfers thermal 
to kinetic energy and leads to vortex formation. However, in astro- 
physical discs, radial stratification is more likely weak. In this case, 
the radial entropy (temperature) variation on the global scale leads 
to the existence of baroclinic perturbations over the barotropic equi- 
librium state. This more appropriate situation has recently become 
a subject of extensive study. 

Klahr and Bodenheimer (2003) pointed out that the radial 
stratification in the disc can lead to the global baroclinic instabil- 
ity. Numerical results show that the resulting state is highly chaotic 
and transports angular momentum outwards. Later Klahr (2004) 
performed a local 2D linear stability analysis of a radially stratified 
flow with constant surface density and showed that baroclinic per- 
turbations can grow transiently during a limited time interval. John- 



son and Gammie (2005) derived analytic solutions for 3D linear 
perturbations in a radially stratified discs in the Boussinesq approx- 
imation. They find that leading and trailing waves are character- 
ized by positive and negative angular momentum flux, respectively. 
Later Johnson and Gammie (2006) performed numerical simula- 
tions, in the local shearing sheet model, to test the radial convec- 
tive stability and the effects of baroclinic perturbations. They found 
no substantial instability due to the radial stratification. This re- 
sult reveals a controversy over the issue of baroclinic instability. 
Presently, it seems that nonlinear baroclinic instability is an un- 
likely development in the local dynamics of sub-Keplerian discs 
with weak radial stratification. 

Potential vorticity production, and the formation and develop- 
ment of vortices in radially stratified discs have been studied by 
Petersen et al. (2007a,b) by using pseudospectral simulations in 
the anelastic approximation. They show that the existence of ther- 
mal perturbations in the radially stratified disc flows leads to the 
formation of vortices. Moreover, stronger vortices appear in discs 
with higher temperature perturbations or in simulations with higher 
Reynolds numbers, and the transport of angular momentum may be 
both outward and inward. 

Keplerian differential rotation in the disc is characterized by a 
strong velocity shear in the radial direction. It is known that shear 
flows are non-normal and exhibit a number of transient phenom- 
ena due to the non-orthogonal nature of the operators (see e.g. Tre- 
fethen et al. 1993). In fact, the studies described above did not take 
into account the possibility of mode coupling and energy transfer 
between different modes due to the shear flow induced mode con- 
version. Mode coupling is inherent to shear flows (cf. Chagelishvili 
et al. 1995) and often, in many respects, defines the role of pertur- 
bation modes in the system dynamics and the further development 
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of nonlinear processes. Thus, a correct understanding of the energy 
exchange channels between different modes in the linear regime is 
vital for a correct understanding of the nonlinear phenomena. 

Indications of the shear induced mode conversion can be 
found in a number of previous studies. Barranco and Marcus 2005 
report that vortices are able to excite inertial gravity waves during 
3D spectral simulations. Brandenburg and Dintrans (2006) have 
studied the linear dynamics of perturbation SFH to analyze non- 
axisymetric stability in the shearing sheet approximation. Tempo- 
ral evolution of the perturbation gain factors reveal a wave nature 
after the radial wavenumber changes sign. Compressible waves 
are present, along with vortical perturbations, in the simulation by 
Johnson & Gammie (2005b) but their origin is not particularly dis- 
cussed. 

In parallel, there are a number of papers that focus on the in- 
vestigation of the shear induced mode coupling phenomena. The 
study of the linear coupling of modes in Keplerian flows has been 
conducted in the local shearing sheet approximation (Tevzadze et 
al. 2003,2008) as well as in 2D global numerical simulations (Bodo 
et al. 2005, hereafter B05). Tevzadze et al. (2003) studied the lin- 
ear dynamics of three-dimensional small scale perturbations (with 
characteristic scales much less then the disc thickness) in vertically 
(stably) stratified Keplerian discs. They show, that vortex and in- 
ternal gravity wave modes are coupled efficiently. B05 performed 
global numerical simulations of the linear dynamics of initially im- 
posed two-dimensional pure vortex mode perturbations in com- 
pressible Keplerian discs with constant background pressure and 
density. The two modes possible in this system are effectively cou- 
pled: vortex mode perturbations are able to generate density-spiral 
waves. The coupling is, however, strongly asymmetric: the cou- 
pling is effective for wave generation by vortices, but not vicev- 
ersa. The resulting dynamical picture points out the importance of 
mode coupling and the necessity of considering compressibility ef- 
fects for processes with characteristic scales of the order or larger 
than the disc thickness. Bodo et al. (2007) extended this work to 
nonlinear amplitudes and found that mode coupling is an efficient 
channel for energy exchange and is not an artifact of the linear anal- 
ysis. B05 is particularly relevant to the present study, since it stud- 
ies the dynamics of mode coupling in 2D unstratified flows and 
is a good starting point for a further extension to radially strati- 
fied flows. Later, Heinemann & Papaloizou (2009a) derived WKBJ 
solutions of the generated waves and performed numerical simula- 
tions of the wave excitation by turbulent fluctuations (Heinemann 
& Papaloizou 2009b). 

In the present paper we study the linear dynamics of pertur- 
bations and analyze shear flow induced mode coupling in the lo- 
cal shearing sheet approximation. We investigate the properties of 
mode coupling using qualitative analysis within the three-mode ap- 
proximation. Within this approximation we tentatively distinguish 
vorticity, entropy and pressure modes. Quantitative results on mode 
conversion are derived numerically. It seems that a weak radial 
stratifications, while being a weak factor for the disc stability, still 
provides an additional degree of freedom (an active entropy mode), 
opening new options for velocity shear induced mode conversion, 
that may be important for the system behavior. One of the direct 
result of mode conversion is the possibility of linear generation of 
the vortex mode (i.e., potential vorticity) by compressible pertur- 
bations. We want to stress the possibility of the coupling between 
high and low frequency perturbations, considering that high fre- 
quency oscillations have been often neglected during previous in- 
vestigations in particular for protoplanetary discs. 

Conventionally there are two distinct viewpoints commonly 



employed during the investigation of hydrodynamic astrophysical 
discs. In one case (self gravitating galactic discs) the emphasis is 
placed on the investigation of the dynamics of spiral-density waves 
and vortices, although normally present in numerical simulations, 
are thought to play a minor role in the overall dynamics. In the other 
case (non-self-gravitating hydrodynamic discs) the focus is on the 
potential vorticity perturbations and density-spiral waves are often 
thought to play a minor role. Here, discussing the possible (multi) 
mode couplings, we want to draw attention to the possible flaws 
of these simplified views (see e.g. Mamatsashvili & Chagelishvili 
2007). In many cases, mode coupling makes different perturbation 
to equally participate in the dynamical processes despite of a sig- 
nificant difference in their temporal scales. 

In the next section we present mathematical formalism of our 
study. We describe three mode formalism and give schematic pic- 
ture of the linear mode coupling in the radially sheared and strati- 
fied flow. Numerical analysis of the mode coupling is presented in 
Sec. 3. We evaluate mode coupling efficiencies at different radial 
stratification scales of the equilibrium pressure and entropy. The 
paper is summarized in Sec. 4. 



2 BASIC EQUATIONS 

The governing ideal hydrodynamic equations of a two- 
dimensional, compressible disc flows in polar coordinates 
are: 

9S | 1 9(rSK) | i £(gy» _ Q; 

dt r dr r d<f> 

m +v m + ndv L _Y± = ^L2f.^^h. (2 ) 

dt dr r d(j> r S dr dr 

dv* { V W4, { v»av» { VrV* _ i dp ^ 

dt dr r d<f> r Sr d<p ' 

dt ^ r dr r d<f> ^ ( r dr r d</> ) ' ^ 

where V r and are the flow radial and azimuthal velocities re- 
spectively. P(r,<f>), T,(r,4>) and 7 are respectively the pressure, 
the surface density and the adiabatic index. ip g is the gravitational 
potential of the central mass, in the absence of self-gravitation 
(ip g ~ —1/r). This potential determines the Keplerian angular 
velocity: 

9ip g 2 -3/2 , c , 

-gjr = Q Ke P r , fiifep ~ r 1 ; (5) 

2.1 Equilibrium state 

We consider an axisymmetric (d/d<p = 0), azimuthal (V r 
f 1 and differentially rotating basic flow: = f2(r)r. In the 2D 
radially stratified equilibrium (see Klahr 2004), all variables are 
assumed to follow a simple power law behavior: 

EM = *,(-) , PV = *>(-) , (6) 

where overbars denote equilibrium and So and P are the values of 
the equilibrium surface density and pressure at some fiducial radius 
r = ro. The entropy can be calculated as: 

S = PS^ = - , (7) 

where 
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Ps = Pp - 



(8) 



S is sometimes called potential temperature, while the physical en- 
tropy can be derived as S — Cv log S + So- 

This equilibrium shows a deviation from the Keplerian profile 
due to the radial stratification: 



1 dP(r) _ 
rE(r) dr 

Po Pp ( r \fe-/3p-2 



= -f±BL. ( L\ 



•(9) 



Hence, the described state is sub-Keplerian or super-Keplerian 
when the radial gradient of pressure is negative (Pp > 0) or 
positive (P P < 0), respectively. Although these discs are non- 
Keplerian, they are still rotationally supported, since the deviation 
from the Keplerian profile is small: Af2 2 (r) <C &K ep - 

2.2 Linear perturbations 

We split the physical variables into mean and perturbed parts: 

E(r,0) = E(r) + E'(r,0), (10) 

P{r^) = P{r) + P'{r,<t } ), (11) 

V r (r,4>) = V;(r,<j>), (12) 

V*(r,0) = fi(r)r + V2(r,0). (13) 

In order to remove background trends from the perturbations we 
employ the global radial power law scaling for perturbed quantities: 



( r ) = fl.y Ss E'(r) , 
V ro / 

(r) V(r) 



Vro 



V'W 



(14) 
(15) 
(16) 



After the definitions one can get the following dynamical 
equations for the scaled perturbed variables: 



-/3 s — <5s+£v 



dV r l + S v -p s 

dr r d(j> r 



l+Sl{r)^}V, -mr)V. + 



c 2 ( r \0e+*p-*v d P +c 2 Sp / r \ 
7 V ro / dr Po 7ro V ro / 



(17) 

= 0, 
(18) 



P, + 



2 /3p 



7ro Vro/ 



V* + 2fi(r) + r 



an(r) 

dr 



V r + 



(19) 



J+5y — <5p 



dr r d<f> r 



where c 2 = 7P0/S0 is the squared sound speed at r = ro. 



2.3 Local approximation 

The linear dynamics of perturbations in differentially rotating 
flows can be effectively analyzed in the local co-rotating shear- 
ing sheet frame (e. g., Goldreich & Lynden-Bell 1965; Goldreich 
& Tremaine 1978). This approximation simplifies the mathemat- 
ical description of flows with inhomogeneous velocity. In the ra- 
dially stratified flows the spatial inhomogeneity of the governing 
equations comes not only from the equilibrium velocity, but from 
the pressure, density and entropy profiles as well. In this case we 
first re-scale perturbations in global frame in order to remove back- 
ground trends from linear perturbations, rather then use complete 
form of perturbations to the equilibrium (see Eqs. 14-16). Hence, 
using the re-scaled linear perturbation (P, E, V) we may simplify 
local shearing sheet description as follows. Introduction of a local 
Cartesian co-ordinate system: 



- r 



y = r ( 



Slot) , -, 2- « 1 
ro r 



(21) 



9 d d Id d d _ d 

r ft ^-, (22) 



dx dr ' dy ro d</> ' dt dt 

where fio is the local rotation angular velocity at r = ro, trans- 
forms global differential rotation into a local radial shear flow and 
the two Oort constants define the local shear rate: 



A=-r 



B 



dn(r) 
dr 

an(r) 

dr 



+ 2fi(r) 



= -A-Sl 



(23) 



(24) 



Hence, the equations describing the linear dynamics of perturba- 
tions in local approximation read as follows: 



d_ 2Ax _d_ 1 P 

dt dy J 7P0 



(25) 



dV x | dV y | 1 + 5 V 
dx dy 



PpIi 
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V x 
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d P | 5 P +p P h P p P S 
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+ 2Ax 



dt 



dy 



Vy - 2BVy + C\ 



d P 

dy 7P0 



[dt dy J 7P0 in) 



(26) 

= 0, 
(27) 
(28) 



where S is the entropy perturbation: 



A. (L 

7r Vro 



\ 0s+S P -6 v -l Q p 
) M^ =0 ' 



§-t + ^-k\T + 



(20) 



s = p - ciz . 



(29) 



Now we may adjust the global scaling law of perturbations in order 
to simplify the local shearing sheet description (see Eqs. 25,26): 



1 + 5 V - PpH = : 



(30) 
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s P + fiph = o 



(31) 



Let us introduce spatial Fourier harmonics (SFHs) of pertur- 
bations with time dependent phases: 

( KM) \ / u x (k(t),t) 

Vy(r,t) = I Uy(k(t),t) 

P(-r,t)hP I -ip(k(t),t) 

\ S(r,t)/ 1 P J V a(k(t),t) 



(32) 



exp (ik x (t)x + ik y y) 



with 

k x (t) = k x (0)~2Ak y t. 



(33) 



Using the above expansion and Eqs. (27-30), we obtain a com- 
pact ODE system that governs the local dynamics of SFHs of per- 
turbations: 



kyUy - 



— p - k x (t)u x 

7-Mj — 2Q. Q u y + c 2 k x (t)p — c 2 s k P s = . 
di 

d 2 

—u y - 2Bu x + c 3 k y p = , 
- k s u x = . 



where 
k P — 



•yr -yr 



The potential vorticity: 

W = k x (t)u y - k y u x - 2Bp , 



(34) 
(35) 
(36) 
(37) 

(38) 

(39) 



is a conserved quantity in barotropic flows: W = const, when 
k P = 0. 



2.4 Perturbations at rigid rotation 

The dispersion equation of our system can be obtained in the shear- 
less limit (A = 0, B = —Q). Hence, using Fourier expansion of 
perturbations in time oc exp(ioii), in the shearless limit, we obtain: 

lu 4 - (c 2 fc 2 + 4^o - c 2 s rj) u 2 - c s r)k 2 y = , (40) 

where 

77 EE kpk S = . (41) 

Solutions of the Eq. (40) describe a compressible density-spiral 
mode and a convective mode that involves perturbations of entropy 
and potential vorticity. For weakly stratified discs (ij <C k 2 ), we 
find the frequencies are: 

Q% = c 2 k 2 + 40?, , (42) 



CgTjky 



(43) 



elk* + 

High frequency solutions (Cu 2 ) describe the density-spiral waves 
and will be referred later as the P-modes. Low frequency solutions 
(d>c), instead, describe radial buoyancy mode due to the stratifica- 
tion. In barotropic flows (77 = 0) this mode is degenerated into sta- 
tionary zero frequency vortical solution. Therefore, we may refer to 
it as abaroclinic mode. The mode describes instability when 77 > 0. 



In this case the equilibrium pressure and entropy gradients point in 
the same direction. Klahr (2004) has anticipated such result, al- 
though worked in the constant surface density limit (/fe = 0). The 
same behavior has been obtained for axisymmetric perturbations in 
Johnson and Gammie (2005). For comparison, in our model baro- 
clinic perturbations are intrinsically non-axisymmetric. Hence, our 
result obtained in the rigidly rotating limit shows that the local ex- 
ponential stability of the radial baroclinic mode is defined by the 
Schwarzschild-Ledoux criterion: 



dPdS 



> . 



(44) 



dr dr 

The dynamics of linear modes can be described using the 
modal equations for the eigenfunctions: 

d 2 



dTi 2 



, -2 



*„,c(t) = ; 



(45) 



where and $ c (t) are the eigenfunctions of the pressure 

and convective (baroclinic) modes, respectively. The form of these 
functions can be derived from Eqs. (34-41) in the shearless limit: 



= (cd 2 >c + c 2 V )p(t) - 2Q W(t) - c 2 k P k x s(t) 



(46) 



All physical variables in our system (p, u x ,u y , s) can be expressed 
by the two modal eigenfunctions and their first time derivatives 
(<& Pi c, 3>p jC ). Hence, we can fully derive the perturbation field of a 
specific mode individually by setting the eigenfunction of the other 
mode equal to zero. 

As we will see later, the Keplerian shear leads to the degen- 
eracy of the convective buoyancy mode. In this case only the shear 
modified density-spiral wave mode eigenfunction can be employed 
in the analysis. 



2.5 Perturbations in shear flow: mode coupling 

It is well known that velocity shear introduces non-normality into 
the governing equations that significantly affects the dynamics of 
different perturbations. In this case we benefit from the shearing 
sheet transformation and seek the solutions in the form of the so- 
called Kelvin modes. These originate from the vortical solutions 
derived in seminal paper by Kelvin (1887). In fact, as it was ar- 
gued lately (see e.g., Volponi and Yoshida 2002), the shearing 
sheet transformation leads to some sort of generalized modal ap- 
proach. Shear modes arising in such description differ from linear 
modes with exponential time dependence in many respects. Pri- 
marily, phases of these continuous spectrum shear modes vary in 
time through shearing wavenumber; their amplitudes can be time 
dependent; and most importantly, they can couple in limited time 
intervals. On the other hand, shear modes can be well separated 
asymptotically, where analytic WKBJ solution for the each mode 
can be increasingly accurate. In the following, we will simply refer 
to these shearing sheet solution as "modes". 

The character of shear flow effects significantly depend on the 
value of velocity shear parameter. To estimate the time-scales of the 
processes we compare the characteristic frequencies of the linear 
modes \Q P \, |u> c | and the velocity shear \A\. In order to speak about 
the modification of the linear mode by the velocity shear, the basic 
frequency of the mode should be higher than the one set by shear 
itself: oj 2 > A 2 . Otherwise the modal solution can not be used to 
calculate perturbation dynamics, since perturbations will obey the 
shear induced variations at shorter timescales. 

In quasi-Keplerian differentially rotating discs with weak ra- 
dial stratification: 
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Gil > A 2 and q\ < A 2 , when < 1 . (47) 

In this case the convective mode diverges from its modal behav- 
ior and is strongly affected by the velocity shear: the thermal and 
kinematic parts obey shear driven dynamics individually. There- 
fore, we tentatively distinguish shear driven vorticity (W) and en- 
tropy (S) modes. On the contrary, the high frequency pressure mode 
is only modified by the action of the background shear. Hence, we 
assume the above described three mode (S, W and P) formalism as 
the framework for our further study. 

For the description of the P mode in differential rotation we 
define the function: 

* p (t) = uj' 2 p (t)p(t) - 2$W(t) - c 2 a k P k x {t)s{t) , (48) 
where 

Wp(t) = c 2 k 2 (t) -4Bfio • (49) 

This can be considered as the generalization of the <6p(t) eigen- 
function for the case of the shear flow, by accounting for the tem- 
poral variation of the radial wavenumber. 

In order to analyze the mode coupling in the considered limit, 
we rewrite Eqs. (34-39) as follows: 

2 

(50) 



d d 



(A , 1 d*„ 

I At + J S = X " pl ~df + Xap2 + Xaw ' 



d_ 

dt 
AW 
At 



(51) 



(52) 



where f p and Alo 2 describe the shear flow induced modification to 
the P-mode 



f P = 4A 



MY 



Ato r , 



kx ky 

k 2 uip ' 

(wjj)" (cu 2 v Y kl 
uii uji k 



parameter f s describes the modification to the entropy mode 



fs = c s ri 



k 2 uip 



(53) 
(54) 

(55) 



and x parameters describe the coupling between the different 
modes: 

k 2 

Xpw = 2fl Auj 2 (t) + -iA-^uj 2 



fc 2 3 

X P s = c 2 k P k x [Au; 2 p +4B^^-8AB^ 



Xspl 



ksk x 

k 2 LO% 



Xsp2 = 



ksk x ( (wp) 



k 2 u>% 



+ 2B- 



2ttk s k x ( (uj p ) 



Xu 



k 2 LOp 



-c„kpk- 



+ 2B-2- + 



kyLUp 



P Ky 



(56) 
(57) 
(58) 
(59) 

(60) 
(61) 



Here prime denotes temporal derivative. 

Equations (50-52) describe the linear dynamics of modes and 
their coupling in the considered three mode model. In this limit, 



our interpretation is that the homogeneous parts of the equations 
describe the individual dynamics of modes, while the right hand 
side terms act as a source terms and describe the mode coupling. 
This tentative separation is already fruitful in a qualitative descrip- 
tion of mode coupling. 

Dynamics of the density-spiral wave mode in the differential 
rotation can be described by the homogeneous part of the Eq. (50). 
Homogeneous part of Eq. (51) describes the modifications to the 
entropy dynamics. Inhomogeneous parts of the Eqs. (50-52) reveal 
coupling terms between the three linear modes that originate due to 
the background velocity shear and radial stratification. We analyze 
the mode coupling dynamics numerically, but use the coupling \ 
coefficients for qualitative description. 

The sketch illustration of the mode coupling in the above de- 
scribed three-mode approximation can be seen in Fig.Q] The figure 
reveals a complex picture of the three mode coupling that originates 
by the combined action of velocity shear and radial stratification. 

The temporal variation of the coupling coefficients during the 
swing of the perturbation SFHs from leading to trailing phases is 
shown in Fig. [2] The relative amplitudes of the \pw and x P a pa- 
rameters reveal that potential vorticity is a somewhat more effec- 
tive source of P mode perturbations when compared to the entropy 
mode. On the other hand, it seems that S mode excitation sources 
due to potential vorticity (\sw) can be stronger when compared 
with the P-mode sources (Xspl. Xsp2). 

The effect of the stratification parameters on the mode cou- 
pling is somewhat more apparent. First, we may conclude that the 
excitation of the entropy mode, which depends on the parapeters 
Xspi' Xs P 2 and \aw is generally a stronger process for higher en- 
tropy stratification scales ks (see Eqs. 58-60). Second, we see that 
the generation of the potential vorticity depending on the Xws pa- 
rameter proceeds more effectively at hight pressure stratification 
scale kp. And third, we see profound asymmetry in the three-mode 
coupling: P-mode is not coupled with the W-mode directly. 

A quantitative estimate of the mode excitation parameters can 
be done using a numerical analysis. In this case, the amplitudes of 
the generated W and S modes can be estimated through the values 
of potential vorticity or entropy outside the coupling area. In order 
to quantify the second order P mode dynamics we define its modal 
energy as follows: 



Ep{t) = \^p(t)'\ 2 + up{t) 2 \^p{t)\ 2 



(62) 



This quadratic form is a good approximation to the P mode energy 
in the areas where it obeys adiabatic dynamics: k x (t)/k y S> 1. 

The presented qualitative analysis suggests that perturbations 
of the density-spiral waves can generate entropy perturbations not 
only due to the flow viscosity (not included in our formalism), but 
also kinematically, due to the velocity shear induced mode cou- 
pling. The generated entropy perturbations should further excite 
potential vorticity through baroclinic coupling. Hence, it seems that 
in baroclinic flows, contrary to the barotropic case, P-mode pertur- 
bations are able to generate potential vorticity through a three-mode 
coupling mechanism: P — > S — > W. We believe that traces of the de- 
scribed mode coupling can be also seen in Klahr 2004, where the 
process has not been fully resolved due to the numerical filters used 
to remove higher frequency oscillations. 



3 NUMERICAL RESULTS 

In order to study the mode coupling dynamics in more detail we 
employ numerical solutions of Eqs. (34-37). We impose initial con- 
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A=0: 



P-mode: <P p , £> p 



Buoyancy: <E> C ,Q C 



A 2 2 

A >co c 



P-mode: ¥ p ,¥ p ' 





Figure 1. Mode coupling scheme. In the zero shear limit two second order 
modes P-mode and buoyancy mode with eigenfunctions <3? p and <J>c are 
uncoupled. In the shear flow, when the characteristic time of shearing is 
shorter then the buoyancy mode temporal variation scale (A 2 > Hi 2 ), we 
use three mode formalism. In this limit we consider the coupling of the P, W, 
and S modes. \ parameters describe the strength of the coupling channel. 
Asymmetry of the mode coupling is revealed in the fact that compressible 
oscillations of the pressure mode are not able to directly generate potential 
vorticity, but still do so via interaction with S-mode and farther baroclinic 
ties with W-mode. 




ditions that correspond to the one of the three modes and use a 
standard Runge-Kutta scheme for numerical integration (MATLAB 
ode34 RK implementation). Perturbations corresponding to the in- 
dividual modes at the initial point in time are derived in the Ap- 
pendix A. 



3.1 W-mode: direct coupling with S and P-modes 

In this subsection we consider the dynamics of SFH when only the 
perturbations of potential vorticity are imposed initially. As it is 
known from previous studies (see Chagelishvili et al. 1997, Bodo 
et al. 2005) vorticity perturbations are able to excite acoustic modes 
nonadiabatically in the vicinity of the area where k x (t) = 0. Here 
we observe a similar, but more complex, behavior of mode cou- 
pling. The W-mode is able to generate P and S-modes simulta- 
neously. Fig. [3] shows the evolution of the W-mode perturbations 
in a flow with growing baroclinic perturbations (rj > 0). The re- 
sults show the excitation of both S and P-mode perturbations due 
to mode coupling that occurs in a short period of time in the vicinity 
of t = 10. The following growth of the negative potential vorticity 
is due to the baroclinic coupling of entropy and potential vorticity 
perturbations. 

Fig. [4] shows the evolution of potential vorticity SFH in flows 
with negative rj. After the mode coupling and generation of P and S- 
modes, we observe a decrease of potential vorticity. This represents 
the well known fact that stable stratification (positive Richardson 
number) can play a role of "baroclinic viscosity" on the vorticity 
perturbations. 

Numerical calculations show that the efficiency of the 
mode coupling generally decreases as we increase the azimuthal 
wavenumber k y corresponding to an increase of the density-spiral 
wave frequency: lower frequency waves couple more efficiently. 

To test the effect of background stratification parameters on 
the mode coupling, we calculate the amplitude of the entropy and 
the energy of the P-mode perturbations generated in flows with dif- 
ferent pressure and entropy stratification scales. The amplitudes are 
calculated after a lOfi^ 1 time interval from the change in sign of 
the radial wave-number. In this case, modes are well isolated and 
the energy of the P mode can be well defined. 

Fig.|5]shows the results of these calculations. It seems that the 
mode coupling efficiency is higher with stronger radial gradients. In 
particular, numerical results generally verify our qualitative results 
that the S-mode generation predominantly depends on the entropy 
stratification scale feg. Therefore, P-mode excitation is stronger for 
higher values of rj. 
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Figure 2. The coupling \ parameters vs. the ratio of radial to azimuthal 
wavenumbers k x (t)/k y when latter passes through zero value during the 
time interval At = 40" 1 . Here k y = k P = k s = 0.5-ff -1 . 



3.2 S-mode: direct coupling with W and P-modes 

Fig.[6]shows the evolution of the S-mode SFH in a flow with grow- 
ing baroclinic perturbations. Here we observe two shear flow phe- 
nomena: mode coupling and transient amplification. Entering the 
nonadiabatic area (around t — 10) the entropy SFH is able to gen- 
erate the P-mode, while undergoing transient amplification itself. 
The transient growth of entropy is unsubstantial and the growth 
rate decreases with the growth of k y . The W-mode is instead con- 
stantly coupled to entropy perturbations through baroclinic forces, 
although higher entropy perturbations at later times give an higher 
rate of growth of potential vorticity. The total energy of perturba- 
tions is however dominated, at the end, by the P mode. 

Fig.|7]shows the dependence of the W and P-mode generation 
on the pressure and entropy stratification scales. As expected from 
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Figure 3. Evolution of the W-mode SFH in the flow with k x (t) = 
-30H- 1 , k y = 2H- 1 andequilibr ium with growing baroclinic pertur- 
bations kp = ks = 0.2H" 1 . Mode coupling occurs in the vicinity of 
t = IOCq 1 , where k x (t) = 0. Excitation of the P and S-modes are clearly 
seen in the panels for pressure (P) and entropy (S) perturbations. Perturba- 
tions of the potential vorticity start to grow due to the baroclinic coupling 
with entropy perturbations. 
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Figure 4. Evolution of the W-mode SFH in the flow with k x (t) = 
— 30-ff -1 , k y = 2H —1 and equilibrium with positive rj: kp = 
— 0.2// -1 , k$ = 0.2H~ 1 . Interestingly, SFH dynamics shows the de- 
cay of potential vorticity after the mode coupling and excitation of S- and 
W-modes at t = lOfl^ 1 . The latter fact is normally anticipated process in 
the flows that are baroclinically stable. 
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Figure 5. Surface graph of the generated S and P-mode amplitudes at 
ky = 2H , kx = — 60-ff - \ and different values of kp and kg. Ini- 
tial perturbations are normalized to set E(0)=1. Excitation amplitudes of 
the entropy perturbations show stronger dependence of the kg (left panel), 
while both entropy and pressure scales are important (approximately kgkp 
dependence) for the generation of P-modes (right panel). See electronic edi- 
tion of the journal for color images. 



Figure 6. Evolution of the S-mode SFH in the flow with k x (t) = 
-30H- 1 , k y = 2H- 1 and equilibr ium with growing baroclinic pertur- 
bations kp = kg = 0.2_ff _1 . Perturbations of the potential vorticity are 
coupled grow from the beginning due to the baroclinic coupling with en- 
tropy perturbations. Excitation of the P-mode is clearly seen in the panel 
for pressure (P), while the panel for entropy perturbations (S) shows swing 
amplification in the nonadiabatic area around k x (t) = 0. Change of the 
amplitude of the entropy SFH affects the growth factor of potential vortic- 
ity SFH. 
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Figure 7. Surface graph of the generated W and P-mode amplitudes at 
ky = 2H~ 1 , kx = —GOH -1 , and different values of kp and kg. Ini- 
tial perturbations are normalized to set E(0)=1. Excitation amplitudes of 
the entropy perturbations show predominant dependence on the kg (left 
panel), while only pressure stratification scale kp is important for the gen- 
eration of P-modes (right panel). See electronic edition of the journal for 
color images. 

qualitative estimates, P-mode excitation depends almost solely on 
the pressure stratification scale kp, while the generation of poten- 
tial vorticity generally grows with rj. 

3.3 P-mode: direct coupling with S-mode and indirect 
coupling with W-mode 

Fig. [8] shows the evolution of an initially imposed P-mode SFH 
in a flow with growing baroclinic perturbations. The oscillating 
behavior of the entropy perturbation for t < 10 is given by the 
P-mode. This oscillating component has a zero mean value when 
averaged over time-scales longer than the wave period. The exis- 
tence of the aperiodic S-mode is instead characterized by a nonzero 
mean value. When the azimuthal wavenumber k y (t) changes sign 
at t — 10, we can observe the appearance of a nonzero mean value 
(marked on the plot by the horizontal dashed line), indicating that 
the high frequency oscillations of the P-mode are able to gener- 
ate the aperiodic perturbations of the S-mode. The aperiodic part 
of the entropy perturbation is than able to generate potential vor- 
ticity perturbations. However, as we see from Eq. (54) and Fig.Q~| 
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Figure 8. Evolution of the P-mode SFH in the flow with k x (t) = 
-30H" 1 , k y = 2H- 1 andequilibr ium with growing baroclinic pertur- 
bations kp = kg = Mode coupling occurs in the vicinity of 
t = 10f2 ( y , where W and S-modes are excited. The amplitude of the gen- 
erated aperiodic contribution to the entropy perturbation is marked by the 
red doted line. Farther, this component leads to the baroclinic production of 
potential vorticity with negative sign. 



there is no direct coupling between P and W-modes. Therefore, the 
P-mode generates the S-mode by shear flow induced mode conver- 
sion, while the W-mode is further generated because of its baro- 
clinic ties with the entropy SFH. We describe this situation as the 
three-mode coupling or in other words, indirect coupling of the P to 
the W-mode. Note, that although the S and W-mode generation is 
apparent from the dynamics of entropy and potential vorticity SFH, 
energetically it plays a minor role as compared to the compressible 
energy carried by the P-mode. 

Fig [9] shows that P-mode generates potential vorticity with a 
positive sign. However, the sign of the generated potential vorticity 
depends on the initial phase of the P-mode. Hence, our numerical 
results show generation of the W-mode with both positive and neg- 
ative signs. 

It is interesting also to look at the P-mode dynamics in flows 
stable to baroclinic perturbations (see Fig.[9). The initially imposed 
P-mode is able to generate the S-mode and consequently the W- 
mode, that gives a growth of the potential vorticity with time. Apart 
from the intrinsic limitations (the dependence of the sign of the 
generated potential vorticity on the initial phase of the P-mode and 
the low efficiency of the W-mode generation), this process demon- 
strates the fact that potential vorticity can be actually generated in 
flows with positive radial buoyancy (77 < 0) and Richardson num- 
ber. 

Fig.UOIshows the dependence of the S and W-mode generation 
on the pressure and entropy stratification scales. In good agreement 
with qualitative estimates, the S-mode excitation depends strongly 
on the entropy stratification scale ks, while the generation of the 
potential vorticity generally grows with r\. 



4 CONCLUSION AND DISCUSSION 

We have studied the dynamics of linear perturbations in a 2D, ra- 
dially stratified, compressible, differentially rotating flow with dif- 
ferent radial density, pressure and entropy gradients. We employed 
global radial scaling of linear perturbations and removed the alge- 
braic modulation due to the background stratification. We derived 
a local dispersion equation for nonaxisymmetric perturbations and 
the corresponding eigenfunctions in the zero shear limit. We show 



Figure 9. Same as in previous figure but for kp = — 0.2H~ 1 and kg = 
0.2H~ 1 . Perturbations are stable to baroclinic forces. However, production 
of the potential vorticity with positive sign is still observed. 
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Figure 10. Surface graph of the generated S and W-mode amplitudes at 
ky = 2H~ l , kx = —dOH^ 1 , and different values of kp and kg. Ini- 
tial perturbations are normalized to set E(0)=1. Excitation amplitudes of 
the entropy perturbations mainly depend on the ks (left panel), while both 
pressure and entropy stratification scales are important for the generation 
of W-mode perturbations (right panel). See electronic edition of the journal 
for color images. 



that the local stability of baroclinic perturbations in the barotropic 
equilibrium state is defined by the Schwarzschild-Ledoux criterion. 

We study the shear flow induced linear coupling and the re- 
lated possibility of the energy transfer between the different modes 
of perturbations using qualitative and a more detailed numerical 
analysis. We employ a three-mode formalism and describe the be- 
havior of S W and P-modes under the action of the baroclinic and 
velocity shear forces in local approximation. 

We find that the system exhibits an asymmetric coupling pat- 
tern with five energy exchange channels between three different 
modes. The W-mode is coupled to S and P-modes: perturbations 
of the potential vorticity are able to excite entropy and compress- 
ible modes. The amplitude of the generated S-mode grows with 
the increase of entropy stratification scale of the background (ks) 
while the amplitude of the generated P-mode perturbations grows 
with the increase of background baroclinic index (rj). The S-mode 
is coupled to the W and P-modes: the amplitude of the generated 
P-mode perturbations grows with increase of the background pres- 
sure stratification scale (kp), while the amplitude of the W-mode 
grows with the increase of baroclinic index. The P-mode is cou- 
pled to the S-mode: the amplitude of generated entropy perturba- 
tions grows with the increase of the background entropy stratifica- 
tion scale. On the other hand, there is no direct energy exchange 
channel from P to W mode and, therefore, no direct conversion is 
possible. Our results, however, show that the P-mode is still able 
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to generate W-mode through indirect three mode P-S-W coupling 
scheme. This linear inviscid mechanism indicates that compressible 
perturbations are able to generate potential vorticity via aperiodic 
entropy perturbations. 

The dynamics of radially stratified discs have been already 
studied by both, linear shearing sheet formalism and direct numer- 
ical simulations. However, previous studies focus on the baroclinic 
stability and vortex production by entropy perturbations, neglecting 
the coupling with higher frequency density waves. 

The most vivid signature of density wave excitation in radially 
stratified disc flows can be seen in Klahr (2004). The numerical 
results presented on the linear dynamics of perturbation SFH show 
high frequency oscillations after the radial wavenumber changes 
sign. However, focusing on the energy dynamics, the author filters 
out high frequency oscillations from the analysis. 

The purpose of numerical simulations by Johnson and Gam- 
mie (2006) was the investigation of the velocity shear effects on the 
radial convective stability and the possibility of the development 
of baroclinic instability. Therefore, no significant amount of com- 
pressible perturbations is present initially, and it is hard to judge 
if high frequency oscillations appear later in simulations. Petersen 
et al. (2007a), (2007b) employed the anelastic approximation that 
does not resolve the coupling of potential vorticity and entropy 
with density waves. Moreover, if produced, high frequency den- 
sity waves soon develop into spiral shocks (see e.g., Bodo et al 
2007). The anelastic gas approximation does intentionally neglect 
this complication and simplifies the description down to low fre- 
quency dynamics. 

Numerical simulations of hydrodynamic turbulence in unstrat- 
ified disc flows showed that the dominant part of turbulent energy 
is accumulated into the high frequency compressional waves (see, 
e.g., Shen et al. 2006). On the other hand, it is vortices that are 
thought to play a key role in hydrodynamic turbulence in accretion 
discs, as well as planet formation in protoplanetary disc dynamics. 
Therefore, any link and possible energy exchange between high fre- 
quency compressible oscillations and aperiodic vortices can be an 
important factor in the above described astrophysical situations. 

Based on the present findings we speculate that density waves 
can participate in the process of the development of regular vortical 
structures in discs with negative radial entropy gradients. Numer- 
ical simulations have shown that thermal (entropy) perturbations 
can generate vortices in baroclinic disc flows (see e.g., Petersen et 
al. 2007a, 2007b). Hence, vortex development through this mech- 
anism depends on the existence of initial regular entropy perturba- 
tions, i.e., thermal plumes, in differentially rotating baroclinic disc 
flows. 

It seems that compressional waves with linear amplitudes can 
heat the flow through two different channels: viscous dissipation 
and shear flow induced mode conversion. However, there is a strict 
difference between the entropy production by the kinematic shear 
mechanism and viscous dissipation. In the latter case, compres- 
sional waves first need to be tightly stretched down to the dissi- 
pation length-scales by the background differential rotation to be 
subject of viscous dumping. As a result, the entropy produced by 
viscous dissipation of compressional waves takes a shape of nar- 
row stretched lines. This thermal perturbations can baroclinically 
produce potential vorticity of similar configuration. However, this 
is clearly not an optimal form of potential vorticity that can lead to 
the development of the long-lived vortical structures. On the con- 
trary, entropy perturbations produced through the mode conversion 
channel can have a form of a localized thermal plumes. These can 
be very similar to those used in numerical simulations by Petersen 



et al. (2007a), (2007b). In this case compressional waves can even- 
tually lead to the development of persistent vortical structures of 
different polarity. Hence, high frequency oscillations of the P mode 
can participate in the generation of anticyclonic vortices that further 
accelerate dust trapping and planetesimal formation in protoplane- 
tary discs with equilibrium entropy decreasing radially outwards. 

Using the local linear approximation we have shown the pos- 
sibility of the potential vorticity generation in flows with both, pos- 
itive and negative radial entropy gradients (Richardson numbers). 
In fact, the standard alpha description of the accretion discs implies 
positive radial stratification of entropy and hence, weak baroclinic 
decay of existing vortices. In this case there will be a competition 
between the "baroclinic viscosity" and potential vorticity genera- 
tion due to mode conversion. Hence, it is not strictly overruled that 
a significant amount of compressional perturbations can lead to the 
development of anticyclonic vortices even in flows with positive en- 
tropy gradients. In this case, radial stratification opens an additional 
degree of freedom for velocity shear induced mode conversion to 
operate. Although, the viability of this scenario needs further inves- 
tigation. 

This paper presents the results obtained within the linear 
shearing sheet approximation. At nonlinear amplitudes, the P mode 
leads to the development of shock waves. These shocks induce lo- 
cal heating in the flow. Therefore, a realistic picture of entropy pro- 
duction and vortex development in radially stratified discs with sig- 
nificant amount of compressible perturbations needs to be analyzed 
by direct numerical simulations. 
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APPENDIX A: INITIAL CONDITIONS 

Here we present the approximations used to derive the analytic 
form of the initial conditions corresponding to individual modes 
in radially stratified shear flows. These conditions are used to con- 
struct the initial values of perturbations in the numerical integra- 
tion of the ODEs governing the linear dynamics of perturbations 
in these flows. We employ different methods for high and low fre- 
quency modes. 



Al P-mode 

P-mode perturbations are intrinsically high frequency and well sep- 
arated from low frequency modes everywhere outside the coupling 
region k x /k y < 1. In order to construct P-mode perturbations we 
use convective eigenfunction derived in the shearless limit and ac- 
count for shear flow effects only in the adiabatic limit: 

tf c (t) = (uj 2 (t) + c 2 V )P(t) - 2n W(t) - c 2 k P k x (t)s(t) , (Al) 

where 

2 t.2 
C 3 T}Ky 



<4(t) = 



c 2 k 2 {t) - 4Sfi 



(A2) 



Although this form of the eigenfunction is not valid function for 
describing W and S modes individually in a sheared medium, it has 
proved to be a good tool for excluding both modes from the initial 
spectrum: 



* c (0) = 0. 



(A3) 



Assuming that we are looking for P-mode perturbations with wave- 
numbers satisfying the condition k x (0) jk y S> 1 we may use the 
zero potential vorticity condition: 



W(Q) = . 



(A4) 



Hence, Eqs. (A3,A4) yield the full set of initial conditions for the 
high frequency P-mode SFH of perturbations: 



p(0) = P , u x (0) = U , 
1 



Uy(Q) 



Mo) 



(k y Uo + 2BP ) 



s(0 ) = ^W±£lzZ Po , 

Cs kp kx (0) 



(A5) 
(A6) 

(A7) 



where Prj and Uo are free parameters corresponding to the two P- 
modes in the system. Specific values of these two parameters define 
whether the potential or kinetic part of the wave harmonic is present 
initially. 



A2 Low frequency modes 

In order to derive the initial conditions for the S and W modes in- 
dividually we employ the second order equation for radial velocity 
perturbation that can be derived from Eqs. (34-37): 

j2 



I- 



+ c 2 s k 2 - ABa 



c 2 s kyW+AAc a k y p .(A8) 



dt 2 

^7 I - li«»u [> ■', =c:k,(t)\\ + 2Bcikr». (A9) 



For low frequency perturbations 
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(A10) 

Assuming that cjc(0) <C c 2 k 2 (0) and neglecting the corresponding 
terms in Eqs. (A6-A7) leads to the following algebraic system: 

[c 2 k 2 - 4Bn ] u x = -c 2 k y W + 4Ac 2 k y p . (All) 
[c 2 k 2 - 4Bn ] u y = c 2 k x (t)W + 2Bc 2 k P s . (A12) 

Hence, we can derive the initial conditions for the low frequency 
modes as follows: 
B 



P(0) = 



2Ac 2 3 k 2 + Bu>$(0) 



(2JWo + c 2 k p k x (0)S ) , (A13) 



u»(0) = ^^^y (-c s k y W + 4Ac 2 k y p(0)) , (A14) 
%(0) = (c 2 k x (0)W + 2Bc 2 s k p So) , (A15) 



where 

2/r.x 2/, 2 



w;(0) = ci(K{0) + k' v ) - 4Bfi . 



(A16) 



Eqs. (A11-A14) give the initial values of perturbation SFHs for S- 
mode when 



W = , So # : 
and W-mode when 
Wo , So = . 



(A17) 



(A18) 



